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ABSTRACT. In this paper we construct energy based numerical methods free of ghost forces in three dimen- 
sional lattices arising in crystalline materials. The analysis hinges on establishing a connection of the coupled 
system to conforming finite elements. Key ingredients are: (i) a new representation of discrete derivatives 
related to long range interactions of atoms as volume integrals of gradients of piecewise linear functions over 
bond volumes, and (ii) the construction of an underlying globally continuous function representing the coupled 
modeling method. 

1. Introduction 

In recent years substantial progress has been made in multiscale modeling of materials, see e.g., EHTOl. A 
class of important problems concerns atomistic-to-continuum coupling in crystals, e.g., the quasicontinuum 
method lfl6l and its variants. Since the continuum model fails to provide an accurate prediction in the 
vicinity of defects and other singularities, coupled atomistic/continuum methods have become popular as an 
adaptive modeling approach over the last years, see e.g. references in lTT2l[T4l[T5ll . The main issue that arises 
in these methods is the proper matching of information across scales. In the first attempts in this direction, ad 
hoc coupling of atomistic and continuum energies resulted in numerical artifacts at the interface of atomistic 
and continuum regions, known as ghost forces, e.g., Q. Therefore, the construction of consistent A/C 
couplings (that are free of ghost forces) is crucial in the numerical modeling of crystalline materials. Further, 
since this problem is one of the better identified mathematical problems related to matching of information 
across scales in materials, it might provide useful insight into the study of multi-scale computational methods 
of a more general nature. 

This paper is devoted to the construction of energy based methods free of ghost forces in three dimen- 
sional crystal lattices. The problem of constructing consistent energies in two dimensional lattices was 
resolved recently by Shapeev lfl4l . see also (H. A key idea in [14] is to express differences (discrete deriva- 
tives) related to long range interactions of atoms as appropriate line integrals over bonds. In two space 
dimensions it is then possible to transform the assembly of line integrals over all possible interactions into 
an area integral, through a counting argument known as the bond density lemma, [14]. This lemma fails 
to hold in three space dimensions, thus the construction of energy based consistent couplings based on 
this approach does not seem to be readily extendable to this case; see [ 15 ] where an interesting attempt to 
circumvent this problem is made. Other papers dealing with similar problems include, e.g., El[T7ll6ll9l[T3l. 

Our work adopts a different approach, based on control volumes associated with bonds, which we call 
bond volumes, and on the construction of an underlying globally continuous function representing the cou- 
pled modeling method. The three dimensional coupled energies constructed in this way are free of ghost 
forces. Moreover, they can be combined in a consistent way to high-order finite element discretizations of 
the continuum region. 
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The paper is organized as follows. In Section 1 we introduce necessary notation. In Section 2 we intro- 
duce suitable finite element spaces and atomistic Cauchy-Born models which are used in the construction 
of the coupled methods. In Section 3 we state and prove a key result, Lemma 3.1, which establishes a 
connection between long range differences and volume integrals of piecewise linear functions defined over 
appropriate decompositions of bond volumes into tetrahedra. In Section 4 we present a conforming coupling 
method based on bond volumes. We note that in the continuum region we use the atomistic Cauchy-Born 
models, introduced in [12]. In Section 5 we show that it is possible to introduce discontinuities at the 
interface, thus allowing greater flexibility in the design of underlying meshes, while still obtaining a con- 
sistent ghost-force-free method. The analysis in this section may lead to the design of more general atom- 
istic/continuum coupled methods based on discontinuous finite elements. Finally, in Section 6 we show that 
one can use finite elements of high order to discretize the continuum region. All methods presented here are 
free of ghost forces; they provide a framework that facilitates the design of several alternative formulations. 

1.1. Notation. Lattice, discrete domain, continuum domain. We let ej be the standard basis vectors for 
R 3 , and choose Z 3 as the three-dimensional lattice. The extension to lattices generated by any three linearly 
independent vectors of R 3 is straightforward since it merely involves compositions with a fixed affine map. 
The scaled lattice is el? = {xi = [xn x ,x^,xi^) = el, £ G Z 3 }, wth lattice distance e = 1/k, k G Z+. 
We will consider discrete periodic functions on Z 3 defined over a 'periodic domain' £. More precisely, let 
Mi G Z+, i = 1,2,3 and define 

:= [-Mi + 1, Mi] x [-M 2 + 1, M 2 ] x [-M 3 + 1, M 3 ]. 
Odiscr := eZ 3 DO, £ := Z 3 n -U. 

£ 

Here Q is the continuum domain; the actual configuration of the atoms is Odiscr, the set of atoms of the 
scaled lattice contained in fl In particular, the convex hull of Odiscr is O. Also £ is the basic lattice period 
in the unsealed lattice Z 3 . 

Functions and spaces. The atomistic deformations are denoted 

ye = y(xi) , x# = el, £ G £ where 

ye = Fx£ + V£, with vg = v{x^) periodic with respect to £. 

Here F is a constant 3x3 matrix with det F > 0. The corresponding spaces for y and v are denoted by X 
and and are defined as follows: 

X:={y:£^R 3 , y e = Fx e + v e , veY, £ G £} , 

V := {u : £ — > M 3 , ue = u(xi) periodic with zero average with respect to £}. 
For functions y, v : £ — > M 3 we define the inner product 

(y,v) e := e 3 ^ y t ■ v e . 

For a positive real number s and 1 < p < oo we denote by W S,P (Q, M 3 ) the usual Sobolev space of 
functions y : Q — > M 3 . By W^ P (Q, M 3 ) we denote the corresponding Sobolev space of periodic functions 
with basic period fi. By ( •, • ) we denote the standard L 2 (Q) inner product; for a given nonlinear operator 
A, we shall denote as well by ( DA, v ) the action of its derivative DA as a linear operator applied to v. The 
space corresponding to X in which the minimizers of the continuum problem are sought is 

X := {y : O — > R 3 , y(x) = Fx + v(x), v G V}, where 

V := {u-.n^R 3 , u G W k ' p {Q.,R 3 ) n W]f{Q.,R 3 ), fudx = 0}. 
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Difference quotients and derivatives. The following notation will be used throughout: 

(1.1) D v y r .= y^- y \ e,£ + r,e£, 

denotes the difference quotient (discrete derivative) in the direction of the vector rj. Also, 

%0(C) := q£ , C = (Ci, C2, C3), 

(1.2) V f 0(C) := {%<KC)}., 

Atomistic and Cauchy-Born potential. We consider the atomistic energy 

(1.3) $ a (y) --=e 3 J2 E <h(PvVt)> 

ees. rieR 

where R C Z 3 is a given finite set of interaction vectors, and the interatomic potential 4> v (-) may vary with 
the type of bond, i.e., ^ may depend explicitly on 77. Further, </>,,(•) is assumed to be sufficiently smooth. 

For a given field of external forces /:£—>■ M 3 , where = f(xf), the atomistic problem reads as 
follows: 

find a local minimizer y a in X of : 

(1.4) 

If such a minimizer exists, then 

(D$ a (y a ),v) e = (f,v) e , for all v ef, 

where 

( J D<I> a (y), U ) £: =e 3 ^ J] ^(Ajltt)^]* 

(1.5) 

We employ the summation convention for repeated indices. 

The corresponding Cauchy-Born stored energy function is E0, 

W(F) = Wcb(F) := E <M Fr ?)- 

Then, the continuum Cauchy-Born model is stated as follows: 

find a local minimizer y in X of : 

(L6) * CB (y)-(f,y), 

where the external forces / are appropriately related to the discrete external forces and 

<S> CB (y) := [ W C B(Vy(x))dx. 
Ja 

If such a minimizer exists, (and is a diffeomorphism on Q,) then 

(1.7) (D$ CB (y CB ),v) = (f,v) , for all v G V, 
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where 

(D$ CB (y), v)= f S ia (Vy(x)) dx = [ S ta (Vy(x)) d a v\x) dx , v G V. 

Here the stress tensor S is defined, as usual, by 

I dFi a J ice 

The stress tensor and the atomistic potential are related through: 

(1.8) $a(F) = d^ F ) = ^ (Fr?) 77c. 

2. Finite element spaces and atomistic Cauchy-Born models 

In the sequel we introduce the finite element spaces used in the rest of the paper. In addition we introduce 
an intermediate model connecting the continuum and atomistic models. We call this the atomistic Cauchy- 
Born model (A-CB). 

Trilinear finite elements on the lattice. Let V £) q be the linear space of all periodic functions that are 
continuous and piecewise trilinear on Q. More precisely, let 

Tq := {K C fl : K = (x h ,X£ 1+1 ) x (x£ 2 ,X£ 2+1 ) x (xe 3 ,xe 3 +i), x £ = (x h ,X£ 2 ,X£ S ) G 0, discr }, 

Ve,Q '■= {v ■ & — > H^ 3 ) v G C(fi) , v\k G Qi(K) and V£ = v(xe) periodic with respect to £}, 

where Q\(K) denotes the set of all trilinear functions on K. Whenever we wish to emphasize that we work 
on the specific cell K = (x£ 1 , X£ 1+ \) x (x£ 2 , X£ 2+ i) x (x£ 3 , Xf. i+ i) we shall denote it by Kg . The elements 
of V £: q can be expressed in terms of the nodal basis functions *&£ = ^(x) as 

V ( X ) = ^2V£^£ 1 (X 1 )^£ 2 (X 2 )^£ 3 (X 3 ), V£ = V(X£), 

test 

where we have used the fact that ^e(x) can be written as the tensor product of the standard one-dimensional 
piecewise linear hat functions ^£ i (xi). Here ^(x^ ) = 6 e .g.. 

For any connected set O such that 



(2.1) O = U KeSq K, 

Sq being a subset of Tq we denote by V £ ,q(0) the natural restriction of V £t Q on O. 

Linear finite elements on lattice tetrahedra. Let V £j t be the space of continuous periodic functions that are 
piecewise linear on lattice tetrahedra. A crucial observation is that there are more than one ways to subdivide 
a given lattice cell K into lattice tetrahedra. Our analysis is sensitive to the choice of such a subdivision. At 
this point we assume that the lattice tetrahedra in the following definition are all of the same type, i.e., they 
have been obtained via the same type of subdivision of each lattice cell. With this in mind, we define 

7t = {T C VL : T is a tetrahedron whose vertices are lattice vertices of K% , xi G Sldiscr}, 

(2.2) _ 

V £>T ■= {v : Q ->■ R , v G C(Q) , v\ T G Pi(T) and v t = v(x£) periodic with respect to £}, 



where Pi(T) denotes the set of affine functions on T. As above, for any connected set O such that O = 
Ut£S t T, St being a subset of Tr, we denote by V £j t(0) the natural restriction of V £ ^t on O. 
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2.1. Atomistic Cauchy-Born models on cells and tetrahedra. A decomposition of the cell Ki with 
a vertex at X£ into six tetrahedra is called a type A decomposition if the diagonals (xe, X£ +ei+e3 ) and 

\%£+e2 i -E£+ei+e2+e3 

) are edges of the resulting tetrahedra, see Fig. 1. In other words, the main diagonal 
), the three face diagonals starting at X£, the three face diagonals starting at X£ +ei+e2+e3 , 
and the edges of Kg, together comprise the edges of the six tetrahedra. Notice that in each tetrahedron 




Figure 1 . A type A decomposition of the cell Kg into six tetrahedra. 

originating from a type A decomposition of a cell, exactly three edges are edges of the original cell, these 
are depicted with solid, black lines in Fig. 2. To define the atomistic Cauchy-Born model on tetrahedra we 
need to define first discrete gradients at each tetrahedron T. To this end, we assume that all cells are divided 
into tetrahedra from a type A decomposition. Let v 6 V £j t- Define X7v as 

(2-3) {Vv\ T \ := D ea v\, 

where the discrete derivatives D ea v\ on the tetrahedron T are just the difference quotients of v along the 
edges of T with directions e a . These are the edges shared with those of Ki, shown in black solid lines in 
Fig. 2. For example, for the tetrahedron of Fig. 2, t> e3 v\ = D ex vU see ( |1.1) , whereas D e2 v\ = D e2 v\ +ei+ez . 
Notice that the definition of these discrete derivatives can be extended to any smooth function. Then for each 
tetrahedron T it follows that 

/ W CB (Vv)dx = —W C b(Vv). 
Jt 6 

Further, let y be a sufficiently smooth deformation. We define corresponding the atomistic Cauchy-Born 
(A-CB) energy 

(2.4) * a > CB (y)--=jY, E E^fe) = iE E wcb&v). 

T£K e (T)rj£R £e£ T&K e (T) 

Now, for a given field of external forces /:£—>■ M 3 the tetrahedra! A-CB problem reads as follows: 

find a local minimizer y a ' CB in X of : 

* a > aB (y a ' OB )-(f,y a ' aB ) £ . 
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%£+ei+e 2 +e 3 



Figure 2. A typical tetrahedron resulting from the decomposition of the cell K%. The 
three edges shown in solid, black lines are also edges of the original cell . 



If such a minimizer exists, then 

(D$ a ' CB (y a ' CB ),v) £ = (f,v) e , for all v 6 Y. 

This atomistic model is consistent, in the sense that the above is satisfied for homogeneous deformations 

(Uf(x) = Fx, x 6 f2): 

(2.5) (D$ a > CB (y F ), v) = 0, y F (x) = Fx, 

for all v £ V £ ,t- To show that, it suffices to observe 

(D^ CB (y F ),v)= £ ^Y, E E *,(VVFi/)-Vt^ 



(2.6) 



E <M-E j E Vwr ? 



^2 4>' V ( F V)- J Vvr]dx = 0. 



An alternative discrete model defined over cells was introduced in lfl2l . The average discrete derivatives 
were defined, e.g., as 



(2.7) 



D ei V£ = - {D ei V£ + D ei V£ +£2 + D ei V£ +e3 + D ei V e + e 2 +e s } • 



This leads to a discrete gradient Vy in analogy to ( |2.3| >; see [12J for details. The corresponding cell atomistic 
Cauchy-Born energy is then defined by 

$a,CB (y) : =e s J2 E ^ = eS E W CB(Vy). 

fe£ r?6R fe£ 
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The corresponding cell atomistic Cauchy—Born problem is 

find a local minimizer y a ' CB in X of : 

* a ' CB (y a ' CB )-(f,y a > CB ) £ . 

This atomistic model is consistent as well, in the sense that 

(2.8) (D$ a ' CB (y F ), v) = 0, y F (x) = Fx , 

for all v G V £ ,q. As before, this is implied by 

(D$*> CB (y F ),v) =e 3 ^ <l>v(yVFV) ' Vwrj 

(2.9) 



V 4>'JFrj) • / Vvr/dx = 0. 
3^ -/n 



It was shown in lfT2l that this model is both energy- and variationally consistent to second order in e, 
approximating the exact atomistic model as well as the continuum Cauchy-Born model. 

3. Bond volumes and long range differences 

To construct methods that couple the atomistic and continuum descriptions we need to relate long range 
differences and derivatives of functions defined over bond volumes. To fix ideas, let rj G R, and define the 
bond as the line segment bg = {x G M 3 : x = xi + t rn < t < 1} with endpoints xg and xe +rj . The set of 
all bonds B v consists of all b = be for £ G £ (but for r\ fixed). For given £ and r\ G 1? with 771 7727/3 7^ 0, 
the corresponding bond volume B^ v is the interior of the rectangular parallelepiped with edges parallel to 
the standard basis vectors e» and main diagonal b%, see Fig. 3. Next we shall establish a connection between 
long range differences and piecewise linear functions defined over type A decompositions of bond volumes 
B^ri into tetrahedra, which is defined in analogy to type A decompositions of cells Kg. To this end let 
Bi^ n {T) a type A decomposition of the bond volume Bi jTI into six tetrahedra, i.e., the decomposition were 
the diagonals (xe, xg +eim+e3V3 ) and (xe +e2ri2 , xe +rj ) are edges of the resulting tetrahedra, see Fig. 3. 

The following lemma plays a central role in our work. 
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Lemma 3.1. Let v be a piecewise linear and continuous function on a type A decomposition of the bond 
volume Bg^ „ into tetrahedra. Then 

(3.1) e 3 D v Vi = -. r / Vv(x)r]dx. 

\VlV2V3\ Jb 1>71 

Proof. We present the proof for t]% > 0, i = 1, 2, 3. The other cases are similar. We have, 
/ Vv(x)r]dx= / vu-r/ds 

V1V2V3 Js Un mmv3 JdB t ^ 

= [ (-r]i)vds+ [ rjivdsX, 

VIV2V3 ^ U«Bi,,(-ei) JdB lin {ei) J 

where dBi^ie-i) is the face of Bg jT) with outward unit normal e^. Since v is linear in each tetrahedron of 
the decomposition of B( „, it will be linear in each of the two triangles comprising the face dB£ v (r]i). 
Therefore, if r is such a triangle, the integral of v over r can be found explicitly: 

f \t\ 3 

(3-3) / rnvds = — V Viv(zj), 

6 j=i 

2 

where Z{ are the vertices of r. Since r is one of the two triangles of dBi^irn), |r| i]i = y i]i i]2 Hence, 

If e 2 2 

(3.4) / wvds = -V {v{z J )+2v(z j )}, 

where Zj are the vertices shared by two triangles of dBg riipi) an d z j tne vertices belonging to only one 
triangle of dBi ^iei). 



We substitute the above formula into ( |3.2| > and group together all terms involving each vertex. For each 
of the vertices other than X£ or X£ +r) , there are two possibilities: 

(i) It is a shared vertex in one face with outward normal e-i and it is a single vertex in two faces with 
normal — e%. 

(ii) It is a shared vertex in one face with normal — a and a single vertex in two faces with normal ej. 
Also, terms involving a vertex of dBi :V (ei) appear with coefficient 1, while terms involving a vertex of 
dBi^ v (—ei) appear with coefficient —1 in ( |3.2| ). Therefore the contribution of these vertices to the sum in 
( |3.2| ) is zero. 

Finally, we notice that xe +v is a shared vertex at each dBi^iei), while xe is a shared vertex at each 
dBi : - v (—ei), for all i = 1, 2, 3. It follows that 

(3.5) / Vv(x) -rjdx = e 2 (v e+v - vA , 

and the proof is complete. □ 

4. A COUPLING METHOD BASED ON BOND VOLUMES 

In this section we construct methods based on bond volumes. Let the atomistic region Q a and the A-CB 
region each be the interior of the closure of a union of lattice tetrahedra T G It and connected, and 
suppose 

si = n a un*, r = n a nn„. 
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Here F is the interface. To avoid technicalities that may arise due to the fact that we work with periodic 
functions over Q,, we assume throughout that Q a is subset of the interior of Q with sufficient distance from 
dtt. Let yi be the deformed position of xe G f^discr- 

Fix 77 G R, with r/i7/2f?3 7^ 0. The cases of degenerate rj can be treated with two and one dimensional 
techniques. We shall construct an energy based coupling method whose design relies on an appropriate 
handling of bond volumes Bi^. We consider three cases depending on the location of each bond volume 

(a) The closure of the bond volume is contained in the atomistic region: v C Q a - 

(b) The bond volume is contained in the region f2* : , cfi». 

(c) We denote by By the set of bond volumes that do not satisfy (a) or (b). In fact, Be, v G By if the 
bond volume intersects the interface: B^ ^ n F ^ or if Bg „ C Q a and Be tV n F 7^ 0. 

If a bond volume intersects dQ, then it is supposed to belong to Q,* by periodic extension. For a fixed r], the 
contribution to the energy corresponding to the atomistic region (case (a)) is 

(4-1) E^ v {y} = e 3 £ ^(D v y e ) . 

_ ees. 

The contribution to the energy from the A-CB region (case (b)) is (cf. ( |2.4[ )), 

(4.2) E^ n {y} = ~ Y, &v(Vyv)= [ <h(yy{x)ri)dx , 

£e£, T£-K e (T), Ten, n * 



y being the interpolant of {ye} in V £i t(^*), see below ( |2.2[ ). 

For each bond volume intersecting F we denote by y e ' v a piecewise polynomial function on B^ v satisfy- 
ing 

i) y^eC(B i>v ). 

ii) Let T(Be :V ) be a decomposition of Bg^ with the following properties: a) if T G T(Bg^) and 
Tc!i» then T is a tetrahedron resulting from a type A decomposition of an atomistic cell K C £7*. 
b) If T G T(Bi jV ) and T C fi a , then T is a lattice tetrahedron. 
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iii) In case ii.b) above, if T has a face on d (B^ ^ n Q a ) \F, then it is part of a conforming decomposition 
that is compatible with decompositions of other bond volumes sharing a face with Bi „. If such an 
attached bond volume is included in Q a , then it is assumed to be type- A decomposed into tetrahedra. 

iv) For T £ T(B i>v ), y l * £ Pi(T) and it interpolates {y e } at the vertices of T. 

Then the energy due to bond volumes intersecting the interface is defined as 
(4.3) £r>}= E W^jL XaM^)dx. 

Remark 4.1. Notice that the energy that corresponds to the bond volume Bi^ „ G By would be 
(4-4) - 1 [ ^(Vy^dx. 

The part of this energy corresponding to Bt tT) n SI* has been already taken into account in E^ cb ^{y} and 
hence it is not included in the definition ofE r ^{y} ■ 

Remark 4.2. The choice of the decomposition 7~(Bg >ri ) and of the associated piecewise polynomial function 
y i,v is somewhat flexible; see Bill for a more detailed discussion. It might even allow vertices that are not 
lattice points. The only essential requirement is that each function defined through T{Bi^ in the proof 



of Proposition 4.1 below, should satisfy v'-" 1 ' 6 H 1 ^) . Depending on the complexity of the interface F one 
can construct such decompositions more or less efficiently. In many cases this can simplify the computation 
of the associated energy E r v {y}- See for example, Figure 4 for such a choice of decomposition. 

We then define the total energy as follows 

(4.5) £bv{y} = J2 £ viy} 

where 

(4-6) £ v {y} = E& atV {y} + E a d f {y} + E T ^{y} . 



4.1. Consistency. The energy ( |4.5| > based on bond volumes is ghost-force free, as we prove in the following 
proposition. 



Proposition 4.1. The energy ( |4.5| > is free of ghost forces, in the sense that 

(4.7) (D£ bv {y F ), v) = 0, y F {x) = Fx , 
for all »£f. 

To show this proposition we shall need some more notation. First we fix rj and consider decompositions into 
bond volumes which cover M 3 : 

(4.8) S% v := [B lv : (i) B e>r) n B,- „ = 0, if £ ± j, (ii) M 2 = } , m = 1, . . . ,\ mmm \. 

will be used for counting purposes in the proof; the associated functions introduced below will be de- 
fined on Q. The number of different such coverings is \r}i 772 % | , hence the numbering m = 1, . . . , \r}i 772 773 1 - 
Notice that bond volumes corresponding to different m may overlap, but the elements of a single S^f are 
non-overlapping bond volumes. 
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For a lattice function \ye} construct the functions Vv and v e ' v in analogy with Vy and y e,v in the con- 
struction below ( |4.2[ ). Then for a fixed 77 we have 

{D£r,{y F ),v} = <t>' v (Fri)-{e 3 Yl D v v t + f Vv(x) V dx 

(4.9) * l "> cn ° 

+ E i 1 / Xu^ndx). 



The main idea in the proof of Proposition 4.1 is to rewrite the expression within brackets above in the 
following way: 



e 3 Yl D v v e + [ Vv(x)r,dx + Y, ■ *„ ■ / 



v e,v r] dx 



(4.10) 



B tyV cn a B t<r) &B T 



\m m ml 



\mv2V3, 

' 1 m=l 



L '14 'IO I /> 

~! J n 



where wl" 1 !, m = 1, . . . , \rji 772 773 1 are appropriate conforming functions in H 1 ^), each associated to a 
different covering Sg consisting of bond volumes. The details are provided below. 



Proof of Proposition 4.1 We use Lemma 3.1 



to write 



Be iV c£l a Bi v cfl a 

where is a piecewise linear continuous function on a type A decomposition of the bond volume B# ^ 
into tetrahedra. The superscript m indicates the covering Sg to which Bg <rj belongs. In fact, can be 
defined globally as follows: For a given lattice function {ve}, a fixed m and a covering f t m l is equal to 

- the piecewise linear interpolant of {vg} on a type A decomposition of the bond volume Bg )Tj into 
tetrahedra if Bg :V C $7 a , 

- 7/^, for 2?^ „ n T 7^ 0, where the piecewise polynomial v ^ on B^ „ is defined through (i-iv) below 

- the piecewise linear function interpolating {vg} at lattice tetrahedra T C Bg v C $7*. 

It is clear by construction that each G i^ 1 (SI) . Further, each tetrahedron corresponds to exactly one 
atomistic cell K C SI* belonging to |r/i 772773) different bond volumes ^, each one belonging to a different 
covering <S^ . Thus for T C SI* we have 

. j |r?i r/ 2 t? 3 I 

(4.12) Vv(x)7]dx = -. r V / Vv [m] (x)rjdx. 

Jt ImmVsl ^ JTnB^eS™ 

Therefore, 

(4.13) / Vv(x)r ] dx = - I r V / Vv [m] (x)7]dx. 

Jo, \ViV2V3\ ^ J a. 



12 



CHARALAMBOS MAKRIDAKIS, DIMITRIOS MITSOUDIS AND PHOEBUS ROSAKIS 



By construction of v e ' v and we have 

j~i \mmm\JB fn 



\mmm\ 

^ E E 



B e , v &B r 



\m mm 



m=l Bt^eS^ '-^.i 



Xn„V« H (x)T7da; 



Thus rewriting ( |4.1 1| ) as 
(4.15) 



n |?7i r? 2 T? 3 1 

^ \vimm 1 



m=l jB t,n 



Vi» H (x)rydx, 



we finally obtain 



jri j-t \mmm\JB e , v 



Bp t v <z£l a 



(4.16) 



1 \vimm\ 

\mv2V3 l ^ 



E 

m=l ^,^5- 

Be, v &Br or B^Cfia 
, 1^71 772 ^73 1 „ 

r V / Vv^(x)r]dx. 

\mmm\ £__^ ./n„ 



m=l 



Hence ( |4.10[ ) follows in view of ( |4.12[ ). Therefore the proof of proposition is complete in view of the 
Gauss-Green theorem. □ 



5. The discontinuous bond volume based coupling method 

In this section we show that is is possible to modify energies to allow underlying functions which might 
be discontinuous at the interface. This allows greater flexibility on the construction of the underlying meshes 
and thus the computation of the energy at the interface might become simpler. To retain consistency the in- 
terfacial energies should include terms accounting for the possible discontinuity of the underlying functions. 
There are many alternatives, such as the possibility of adding extra stabilization terms, compare to [ 1 ]. The 
purpose of this paper is however to present the general framework and we will not insist on the various 
modifications and extensions of the methods developed herein. 

Let S7, , il* and V be as in the previous section. Further, we distinguish the same cases a), b) and c) 
regarding the location of each bond volume B^ v . The corresponding energies are still defined by 

(5-1) E^ v {y} = e 3 £ ^(PnVt), 

_ fe£ 
Be v <z£} a 

and 

(5.2) E %*{y} = j Y, <t>v(yyv)= f <j> v (Vy( x )v)dx, 

te£, T£K e (T), Ten, n * 
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y being the piecewise linear function at the lattice tetrahedra interpolating {ye}. The main difference to the 
previous construction in Section 3 is the choice of y i,v and the corresponding energies for each bond volume 
intersecting the interface. In fact we let 

i) ^ec(%\r). 

ii) Further, let T(Bi : v ) be a decomposition of B^ v with the properties a) if T G T(Bi :V ) and Tc !]» 
then T is an atomistic tetrahedron resulting from a type A decomposition of an atomistic cell, b) If 
T £ T(Be jV ) and T d£l a then T is a lattice tetrahedron. 

iii) In the case ii.b) above if T has a face on d(B£ >v n f2 a )\r then it will allow for a compatible 
conforming decomposition with respect to attached bond volumes. In that case if the attached bond 
volume is included in Q a it is assumed to be type-A decomposed into tetrahedra. 

iv) For T £ T(B e>r) ), y*** G Pi(T), interpolating {y e } at the vertices of T. 

We have kept the same properties, but we allow discontinuous matching across the interface T. This provides 
greater flexibility on the construction of y e,v since it allows the presence of arbitrary hanging nodes on the 
interface of the two regions. 

We then define the energy due to bond volumes intersecting the interface as 

(5.3) Ej? tn { y }= r^—\\f XnMVy"' v n)dx- [ </>' v aW' n v})iy e ' v v}dS 

Here, [to?]], {{u>}} denote the jump and the average of a possibly discontinuous function on the interface 

(5.4) {wrjl := (vn a ■ rj) w~ + ■ rj) w + , {w} := ^{uT + w + } , 

where w~ and w + being the limits taken from Q a and respectively, and UQ a , z/^ the coiTesponding 
exterior normal unit vectors, with UQ a = —vn, on T. 

A key observation here is that {y} does not induce inconsistencies on the energy level. In fact, it is 
obvious that if y e,ri G C(Bg„) as in the previous section, then 

(5-5) E° v {y} = E FtV {y} , 

since the extra term on the interface vanishes. Then, as in the previous section, we define the total energy as 
follows: 

(5.6) ££,{y} = Y, £ «{y}> 

where 

(5-7) £°{y} = E& a;ri {y} + Eg v {y} + E° v {y} . 

Despite the fact that we allow discontinuities, the energy £® v is still ghost-force free: 



Proposition 5.1. The energy (5.6 1 is free of ghost forces, in the sense that 
(5.8) (D£? v {y F ),v) = 0, y F (x) = Fx, 

for all v G T. 



Proof. The structure of the proof is the same to that of Proposition 4. 1 hence we present in detail only the 
main differences. We still need the coverings S r g and recall that their elements define a decomposition of 
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non-overlapping bond volumes. As in the proof of Proposition 4. 1 for a given lattice function {vg} we define 
the functions Vv and v ,r> in analogy with Vy and y i ' v , cf., ( |4.2[ ). Then, we have 



(D£°(y F ),v) = ^(Fri) ■ [e 3 £ D v vg + f Vv(x)r, 

Or- O J 



dx 



(5.9) 



Indeed, to show this it suffices to evaluate 



(5.10) 



which is equal to cf>' (Frj) ■ {vrjl for w = y F . 



In parallel to the proof of Proposition 4. 1 we shall prove 



(5.11) 



1 

?/2 f?3| 



f 



m=l 



VvW(x)rjdx- / {v^^jdS 



where i;[ m l , m = 1, . . . , \rji r/2 773 1, are appropriate functions in i/ 1 (il\r)nC(ri\r), possibly discontinuous 
at T; each is associated to a different covering S"g consisting of bond volumes. Relation ( |5.11| > then 
implies (D£® (y F ),v) = since by the Gauss Green theorem, 



Vv [m] (x)r]dx= / {v [m] rjldS. 



It remains therefore to establish ( |5.11| ). To this end, we proceed exactly as in the proof of Proposition 4.1 
In particular, for a given lattice function {ve}, a fixed m and a covering S^ 1 define as 

- the piecewise linear interpolant of {vg} on a type A decomposition of the bond volume Bg^ into 
tetrahedraif Bg v C 

- v e,v , for Bg iV (IT / 0, where the piecewise polynomial on Bg iV , v ,<n is possibly discontinuous on 
Bg^ n T, and is defined through (i-iv) above, 

- the piecewise linear function at the lattice tetrahedra interpolating {vg}, if T is an atomistic tetrahe- 
dron such that T C Bg :V C SI*. 
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Now, by construction G H l (fi\r) n C(0\r), and is possibly discontinuous at F. The rest of the proof 
is identical to the one of Proposition |4. 1 1 with the exception that ( |4.14| ) should be replaced by 



f-i Wimm Jb, „ JB f „nT 



fe£ '** /J| Ja t,v 



(5.12) 

—r E E / B X^(x) v dx- lv H v]dS} 



\VlV2V3\ m=1 B itV ES™ JD i,^ 

Be,r,£B r 



with ( |4.16| > modified accordingly. □ 

6. High-order finite element coupling 

In this section we shall see how the previous analysis can lead to energy-based methods which employ 
high-order (even hp-) finite element approximations of the Cauchy-Born energy on the continuum region 
while remaining ghost-force free. To this end let T ac be a decomposition of $7 into elements with the 
following properties: Let S7, S7 a and S7* as before with F the interface. The approximations will be based on 
decompositions of the continuum region SI* that are compatible on F to V £ ,t{^*)- To this end, let 

7^(S7*) be a conforming decomposition of SI* into tetrahedra with vertices lattice points, 
(6.1) _ 

such that, if T G %(£l*),T (IF + 0, then T G 7^(0*) ■ 

We consider the discrete space 

V M c(0*) = K C(0*) : v\ T 6 Pi(T)forTnr ^ and 

<6 ' _) u| T G Pfc(r) for all other T G 7^*)}. 

This space can be extended to include the atomistic region as well by 

V Kac = {d£ C(H) : u| T S Pi(T) for T G T T (S7 a ) and 

(6.3) 

w|t £ ^h,ac(^*) f or all 7" G 7^(S7*), u periodic on $7 } . 

For v G V^ qc one can define the corresponding lattice function {vf\, simply by interpolating. Conversely, 
for given {vi\ one can find v G V \,ac that coincides with corresponding values of {vg } at the vertices. 
However, at the regions using high-order finite elements the other degrees of freedom should be defined 
with some care. In the following, we assume that we are given a function y, y{x) = Fx + v (x), v G Vh,ac 
and we shall define its atomistic/continuum energy. To this end, 

(6-4) £ h = V £h, v , 



veR 



where 



(6-5) £ h , ri {y} = {y} + / ^{VyixWdx + E F {y} . 

Here for an atomistic point xt G Q a , ye = y(xe), and the local energies Eq „{y}, E T {y} are defined as 



in Section 4, see (4.1 1, (4.3 1. 



The above method can designed to be of arbitrary high order accuracy of the Cauchy-Born energy at 
the continuum region SI* . Such methods are of importance since, by tuning the discretization parameters 
(decomposition of SI* and polynomial degrees) we have the possibility of matching the ideal accuracy at 
the continuum region which is 0(e 2 ). The energy £h is ghost force free. 
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Proposition 6.1. The energy ( |6.4[ ) is free of ghost forces, in the sense that 
(6.6) (D£ h (y F ),v) = 0, y F (x) = Fx, 

for all v G V h>ac . 
Proof. Since v 6 ac we have 



(D£h,r,{yF),v) = 0: 



;(Fr/)-{e 3 ^ ZV^+ / V«(x)r ? dx 



(6.7) 



_ £e£ 

?) CSl a 

+ E 



1 



|f7l*72*73| JB 



Xn^^dx} 



Exactly as in the proof of of Proposition 4. 1 we may write the first and the third term of the the above sum 
as 

1 



e 3 Yl D v v ? + E 



Xn n V?J ,V V dx 



e, 7] 



(6.8) 



\m m m\ 



- | •/ J. • I& '/O I « 

1 r V / Vv [m] (x)r]dx 

mm ml ^ Jn a 



where v^ m \m = 1, . . . , \rji 772 773 1 are the functions defined in the proof of of Proposition 4.1 Define now, 

f I Tp\mmv3\ v [m]( \ for r G Q 

(6.9) S(x) = ^ hi^r,3|^™=l v W» torxGiia, 

[ u(a?), for if!!,. 

Then tracing back the definition of ac at the elements next to the interface T and the proof of Proposition 



4.1 we can show that v is continuous at the interface T. Thus v G H 1 ^), and is periodic. Further, by the 



Gauss-Green theorem, 



(6.10) (D£ Kv (y F ),v) = 4>'{Frj)- { f Vv V +[ V« r?l = fiJFr,) ■ f Vv r, = , 

and the proof is complete. □ 
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